#This file converts only one year of data as an example.

# total so2 emissions

setwd("D:/Replication/TOTALS_nc") ##set up your working directory
rm(list=ls())

library(ncdf4)
library(maps)
library(data.table)
library(foreign)

# open 2015 emissions data
ncin <- nc_open(file.choose(), verbose = TRUE, write = FALSE)

# change name
dname<-"emi_so2"

# get longitude and latitude
lon <- ncvar_get(ncin,"lon")
nlon <- dim(lon)
head(lon)
lat <- ncvar_get(ncin,"lat")
nlat <- dim(lat)
head(lat)
print(c(nlon,nlat))
# get time
nt <- 1

# get emi_so2
emi_so2_array <- ncvar_get(ncin,dname)
dlname <- ncatt_get(ncin,dname,"emi_so2")
fillvalue <- ncatt_get(ncin,dname,"_FillValue")
dim(emi_so2_array)

# reshape the array into vector
emi_so2_vec_long <- as.vector(emi_so2_array)
length(emi_so2_vec_long)
# reshape the vector into a matrix
emi_so2_mat <- matrix(emi_so2_vec_long, nrow=nlon*nlat, ncol=nt)
dim(emi_so2_mat)
head(na.omit(emi_so2_mat))
# create a data frame
lonlat <- as.matrix(expand.grid(lon,lat))
emi_so2_df02 <- data.frame(cbind(lonlat,emi_so2_mat))
names(emi_so2_df02) <- c("lon","lat","emi_so2")
head(na.omit(emi_so2_df02, 20))
head(na.omit(emi_so2_df02))
dim(na.omit(emi_so2_df02))

# restrict data to China
emi_so2_df02$country<-NA

emi_so2_df02$country <- map.where(database="world", emi_so2_df02$lon, emi_so2_df02$lat)
CN <- subset(emi_so2_df02, country=="China")

mydata <- data.frame(CN$lon, CN$lat, CN$emi_so2)

setnames(mydata, old=c("CN.lon","CN.lat","CN.emi_so2"), new=c("lon", "lat","so2"))
write.dta(mydata, "so2_total_2015.dta")
